Statistical properties of fractures in damaged materials. 
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We introduce a model for the dynamics of mud cracking in the limit of of extremely thin layers. 
In this model the growth of fracture proceeds by selecting the part of the material with the smallest 
■ (quenched) breaking threshold. In addition, weakening affects the area of the sample neighbour to 

^\ | the crack. Due to the simplicity of the model, it is possible to derive some analytical results. In 

. particular, we find that the total time to break down the sample grows with the dimension L of the 

t-H ' lattice as L 2 even though the percolating cluster has a non trivial fractal dimension. Furthermore, 

p- j , we obtain a formula for the mean weakening with time of the whole sample. 

a ■ 

Among the phenomena leading to formation of cracks in materials 0], fractures arising from desiccation are 
widespread both in nature where mud cracks of different length scale are present in the landscape J^] and in the 
industry, where paint drying and concrete desiccation are intensively studied ||@]- Despite its relevance, the mecha- 
I ' nism and the dynamics of cracking by desiccation is, however, still unclear. Recently, an experimental work (5) shed 
. some light onto the problem. According to this work, the main source of stress is given by the local friction of the 
C ' layer of material with the bottom of the container. Moreover, the characteristic scale of the crack patterns varies 
linearly with the layer thickness. In the limit of zero thickness the crack patterns loose their polygonal structure (the 
characteristic size of the polygons is zero) and become branched fractals. 

In this letter we introduce an extremely simplified lattice model for cracking of thin mud (or paint) layers, inspired 
Q\ ' by the vectorial and scalar models described in [HJ7) and by invasion percolation [|| , for which a theoretical approach, 
0^ . based on the Run Time Statistics (RTS) scheme [p), allows us to compute some relevant quantities, such as the 
evolution of the breaking probability. Many cellular automata models for crack propagation present as a basic 
ingredient a non-local Laplacian field (electric field, electric current) which drives the formation of the cracks, acting 
on a solid network of bonds or sites (7). In others the stress field evolves by keeping minimum the energy of the 
system. In such a case the components of the vectorial equations obtained are similar to the equations describing the 
£^ ' action of a Laplacian field Q . 

In this model no explicit Laplacian field is present, since the stress is represented by a local random breaking 
J-^ ' threshold. At each time step, the bond with the smallest threshold is removed from the lattice. Short ranged 
correlations are introduced through a damaging of the thresholds of the bonds nearest neighbours of the bond removed. 
In defining this model we are driven by the above cited experimental observation |j] that in an extremely thin layer 
of mud or paint the only source of stress is the local friction with the container. No globally applied Laplacian field 
seems to be present. Moreover, since the drying mud is a mixture of a liquid and a solid (usually amorphous) phase, 
no long range stress relaxation is present, although the growing crack can affect the properties of the medium in its 
vicinity. The simplifications with respect to Refs. restrict the applicability of the model to the special cases above 
mentioned, yet it allows one to write down explicitly some equations for the evolution of the breaking probability. 

The model is defined as follows: On a square lattice a quenched random variable Xi is assigned to each bond i 
where the x^s are extracted from a uniform distribution between and 1. The bond with the lowest value of the 
variable is selected and removed. Then damage is applied, and the bonds nearest neighbours are weakened, i.e. a 
new threshold x\ is extracted with a uniform probability between and the former value Xi. This should mimic the 
stress enhancing nearby crack tips [Q. Then, the next bond with the minimum value of threshold is removed from 
the system and a new weakening occurs. The process continues until a percolating cluster of fractures divides the 
sample in two separate parts. ^From the point of view of mud cracking, the two-dimensional lattice represents a very 
thin layer of mud (paint), and quenched disorder accounts for local stress induced by inhomogeneous desiccation of 
the sample. The dynamics is assumed to be quasi-static, since we assume evolution of cracks in mud desiccation to 
be a slow process. Some authors correctly point out that otherwise time-dependent effects and a non-equilibrium 
dynamics are relevant in crack propagation JToj] . 

In this model we have eliminated the explicit presence of an external field (applied stress) and of the response of 
the material (strain of bonds). The only quantity present is the breaking threshold, whose dynamics is chosen to 
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reproduce the evolution of cracks. This simulates the presence of a local stress field, acting not on the boundaries but 
directly on each bond. Our assumption is based on the experimental results in ref. where as the mud layer becomes 
thinner only the inhomogeneities drive the nucleation of cracks. Furthermore the hypothesis of crack developing under 
the same state of strain not only is usually applied in the presence of thermal gradients [ fi"l| , but is also commonly 
reported in experiments of loading of softened material |J. Hence, such a model is particularly suitable to describe, 
for example, paint drying, where the stress applied to the painted surface depends on the local action of external 
condition (density gradient in the paint). 

Despite the simplicity of the model rule, the findings are rather interesting. We have performed numerical sim- 
ulations, with cylindrical symmetry (periodic boundary conditions in the horizontal direction) for various system 
sizes L. The dynamics stops as soon as a crack spanning the system in the vertical direction appears. In Fig.|l| a 
typical percolating cluster is shown. The fractal dimension of the percolating cluster is computed with a box-counting 
method. We find Dt = 1.77 ± 0.02 for all values of L. We restrict the analysis to the spanning cluster to reduce the 
finite size effects present for the smaller clusters. Also the distribution of finite clusters is not trivial, showing a clear 
power law with exponent r c = 1.54 ± 0.02 (see FigJ^a). 

A recent interest has focused on acoustic emission in order to understand if the power law observed both in 
experiments [[l2| and in models J^J^] is related to Self-Organized Criticality (SOC) fl3|| . The presence of SOC would 
mean that the dynamics of fractures leads the sample to a steady state where small variation of the external field 
can trigger reaction at any length-scale. In particular the external field in this case is the applied stress, and the 
response of the sample can be considered as the energy released (acoustic emission) by one avalanche of cracks, where 
avalanche means a causally connected series of breakdowns. In this over-simplified model the external stress can be 
considered constant, since the only change after any single breakdown is the damaging of the nearest neighbours. 
Then we monitored, as a measure for the acoustic emission the size of an avalanche which is defined as follows. Let 
us suppose that a bond i grows (i.e. it is broken) at time i; this is the initiator of an avalanche, which is defined as 
the set of events geometrically and causally connected to the initial one (bond i). "Causal" connection refers to the 
weakening following any bond breaking. In particular, when bond i grows at time i, the avalanche goes on at time 
t + 1 if a non-broken first neighbour bond j of i is removed. At time t + 2 the avalanche goes on if a bond k grows 
where A: is a non-broken first neighbour of i or j and so on. In FIG.|b we show a linear-log plot of the probability 
distribution of avalanche size, versus sample size L. We get, after a power law transient, an exponential distribution 




FIG. 1. A typical spanning fracture ob- FIG - 2 - ( a ) Probability* distribution (logio-logio 

tained with the model. In the simulations P lot ) A=(s) of th e size of the clusters vs the size, for 

we deal with periodic boundary condition L = 64 > 128 > 256 together with a least square fit. (b) 

on the sides, i.e. left and right boundary Avalanche size distribution (linear-Zogio plot) D(s) 

coincide f° r L = 1 28 > 2 56. While the cluster distribution is 



power law, the avalanche distribution is exponential. 

signalling that a characteristic size exists for the avalanches. This result for avalanches is similar to those obtained 
for a scalar model of Dielectric Breakdown, but differs from the avalanche behaviour in models of fracture . The 
reason for that is twofold. Firstly, in the present definition of an avalanche the threshold is changed only for the 
nearest neighbours. This introduces a typical length scale, while other definitions consider as the threshold the ratio 
between local field and resistivity, thus giving the possibility of large scale correlations. Secondly, in this model grown 
bonds are removed from the system. This represents a substantial difference with the SOC models presented in the 
literature. For example, in a simple toy model of SOC due to Bak and Sneppen |l4| (where a similar refresh for the 
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threshold is present) the dynamics produces clear power laws in the avalanche distribution. There, each site (species) 
deleted is replaced by a new one and is not definitively removed. This is a crucial point, since indeed power law 
behaviour in the presence of a scalar field seems to be related to a "reconstructing rule" that allows one to deal with 
a system where bonds are not removed but they change their status. Then only in the case of plastic deformation, 
one is in the presence of a steady state, as correctly pointed out by Ref. Q. 

Since the evolution of the crack is taken to be extremal, i.e. the bond with the minimal threshold is removed from 
the system, one is allowed to apply a recent theoretical tool, the Run Time Statistics (RTS) to this model in order 
to predict the probability of breakdown of one bond and then the resistance of the whole medium. RTS has been 
introduced to map deterministic quenched extremal dynamics to a stochastic description. As discussed in 0] we can 
associate with every bond j an effective probability density Pj,t(x) for the threshold Xj conditioned by the growth 
history. We can write the growth probability fii, t of the bond i at time t as 
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In this case (contrary to the case of Invasion Percolation), dCt is the whole set of non-broken bonds. Note that at 
time t the number of bonds in dCt is (2L 2 — t), i.e., the total number 2L 2 of bonds in a square lattice with side L 
minus the number of broken bonds before time t. Now we can update every pj, t (x) by conditioning them to the latest 
growth event, this giving the pj t t+i(x)'s. We call mi tt +i(x) the effective probability density at time t+1 of the latest 
grown bond i. It is given by 
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For the remaining bonds we have to distinguish between the weakened ones (the nearest neighbours of i) and the 
rest. For the latter we have 
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Pj,t+i(x) = PjA x ) dyp lyt (y) 
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while for the former we obtain 

Pj,t+i(x) 
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Eqs. |I| |[ ||) coincide with the ones introduced for the RTS approach to Invasion Percolation (apart from the 
different definition of the growth interface dCt). Eq. (||), instead, is new and refers to the nearest neighbours 
weakening. Eqs. ([y g, ||, |^) allow one to study the extremal deterministic dynamics as a stochastic process. In 
particular, \ii yt can be used to evaluate systematically the statistical weight of a fixed growth path, while pj t t{x) stores 
information about the growth history. A very important quantity to characterize the properties of the dynamics is 
the empirical distribution (or histogram) of non-broken thresholds. This quantity is defined as: 



h t (x) 



jedCt 



(5) 



where, h t (x)dx is the number of non-broken bonds between x and x + dx at time t. 
Considering the effect of the growth of bond i at time t on this quantity, 

hj^x (x) = ht{x)- mjj+i (x) pj, t (x)+y^ j pj tt+ i(x) 



(6) 



where j(i) indicates the set of non-broken bonds j nearest neighbours of i and 771^4+1(2;) and pj_ t +\{x) are given 
respectively by Eq.(0) and Eq.(EI). In order to evaluate the statistical properties of the crack evolution, we averaged 
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over all the possible path growths until time t + 1. Introducing the notation (...) for this average, the l.h.s. of Eq.(g) 
can be computed as 

(ht +1 (x)) = \\dC t+1 \\0 t+1 (x) = {N-(t + l)} ( fi t+1 (x), (7) 

where N — 2L 2 is the total number of bonds in the lattice and (fit (x) represents the average thresholds distribution 
over the non-broken bonds at time t (normalized to 1). For the r.h.s. of Eq.(^|) the main difficulty arises in the 

evaluation of (rriij+i) and Pj.t+i( x )^- Following || we can write 

(m itt+1 ) ~ (N-t)<fi t (x)e- {N - t) f° dvMv) . (8) 

Assuming x to be independent of the number n t of bonds weakened at time t and applying the same mean field 
approximations used to obtain Eq.(|8|) P] we have 

( ^2pjA x ) ) = n tM x )' ( 9 ) 

while for the annealed bonds, we obtain, after some algebra, the expression 



Then we can write the following equation for the (fi t +i(x): 



l_e-vf'-wj (io) 



+ n t ^ ¥ jj^Mll-e-^S:^ (11 ) 

Note that even at percolation time N ~ t is a large number. For this reason terms in eq.([ll]) which contain the 
exponentials are negligible for those x for which f Q dy<fi t {y) is finite (larger than l/(N — t)). It is easy to show that 
the continuum limit of Eq.(|ll]) for such values of x is invariant under the rescaling L — » aL (i.e. N — » a 2 N) and 
t —* a 2 t. This result is based on the assumption that: 

n t (L) =n a 2 t {aL) (12) 

The numerical simulations suggest the following scaling form for n t (L): 



n t (L) = n r , 



l + t/AL 2 



(13) 



where (5 = 0.23 ± 0.02 and n max — 6 is the lattice coordination number. This form for n t satisfies Eq. (|l2|) . 

We are then able to make three important predictions: Firstly, we find both analytically, from the numerical 
solution of Eq.jll], and from computer simulations a discontinuity in the histogram (see Fig.^J) , signalling that the 
i system evolves in a way such as to remove all bonds with threshold smaller than some critical value. Secondly, 
from the symmetry properties of Eq. [ll] we deduce that the number t sp (L) of broken bonds at the percolation time 
is proportional to L 2 , even though the percolating cluster is fractal. This result is deduced supposing that at the 
percolation time the shape of the histogram is independent of L which is verified by numerical simulations (Fig.|^-a) . 
The numerical check of this result is presented in Fig.^a. Thirdly, we present a qualitative result for the average value 
of the thresholds (i.e. the resistance of the material) (x) (t). For discretized times we find: 

<*>(* + 1) - (i - 2( ;r t 2 _ 1} ) <*>(*) + 1+n ^ N ; t)] x Jl dxe-^i:^\ ^ 

It is possible to show that {x)(t + 1) < (x)(t) unless (x)(t) is exponentially small in N — t. TOhis means that on 
average the medium weakens during the evolution even if the weakest bond is removed at any time step. In Fig^b we 
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■ ■ num. simulation, L=32 

*■ - - « num. simulation, L-64 
• — * num. simulation, L— 128 
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FIG. 3. Solution of EqJTl] for the histogram 
4>t(x) at the spanning time (b), compared with 
simulations (a). 



■ numerical simulations L-128 
theory, L=128 
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FIG. 4. (a) Spanning time versus system size 
L. One can see a nice agreement with the ex- 
pected scaling law t sp (L) oc L 2 . (b) Solution of 
eq.[14| compared with numerical simulations. 



show the time evolution of (x)(t) obtained from computer simulations, compared with the solution of eq. (J14]) . Our 
analytical results are in nice agreement with numerical simulations. 

As a conclusion, we would like to point out that, to our knowledge, apart the qualitative results of [||, no quantitative 
experimental results are available. For example, a measurement of the fractal dimension of cracks or their size 
distribution would be extremely useful to further test the predictions of this model. At the moment, this model seems 
able to capture, with its extremely simplified dynamics, some basic properties of fracturing processes. 

In conclusion we have presented a new model for fractures, that we believe should be applicable to the case of 
drying paint and mud cracking, for extremely thin samples. Due to its extreme simplicity, the model is particularly 
suitable for large scale simulations and takes into account the damaging effects involved in fracture propagation. Even 
in this simple model we are able to analyze which conditions trigger SOC behaviour in such systems. Furthermore, 
the change in the threshold distribution induced by the annealing allows us to write down explicitly the form of the 
breakdown probability for the bonds of the sample. In the future we plan to extend the analytical study of the model 
to the computation of its the relevant critical exponents. This is in principle possible by combining the RST method 
with the fixed scale transformation approach (FST) fll5|j . The authors acknowledge the support of the EU grant 
Contract No. FMRXCT980183. GC acknowledges the support of EPSRC. 
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